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Parallel  channel  configurations,  such  as  Z-type,  used  to  distribute  reagents  in  planar  fuel  cells  provide 
lower  overall  pressure  drop  as  compared  to  other  channel  designs.  However,  due  to  their  inherent 
characteristics,  flow  maldistribution  in  parallel  configurations  is  commonly  observed  and  leads  to  star¬ 
vation  of  reagents  in  middle  channels.  In  addition,  the  Reynolds  number  dependent  minor  losses  at 
branching  tee  junctions  may  cause  asymmetric  flow  non-uniformity  and  reagent  imbalance  between  the 
cathode  and  anode.  Herein,  we  present  a  universal  and  simple  optimization  method  to  simultaneously 
reduce  flow  maldistribution,  asymmetry,  and  parasitic  pressure  in  Z-type  parallel  configurations  of  fuel 
cells  or  fuel  cell  stacks  that  has  improved  scalability  relative  to  previous  methods.  A  discrete  model's 
governing  equations  were  reduced  to  yield  geometric  ratios  between  headers.  Increasing  header  widths 
to  satisfy  these  ratios  reduced  flow  maldistribution  without  modifying  parallel  channel  geometry  as 
validated  by  computation  fluid  dynamics  (CFD)  simulations.  Furthermore,  decreased  Reynolds  numbers 
throughout  the  headers  reduced  minor  pressure  drops  and  flow  distribution  asymmetry.  We  offer 
several  methods  to  reduce  the  optimized  geometry's  footprint,  including  an  adaptation  of  the  discon¬ 
tinuous  design. 

©  2014  Elsevier  B.V.  All  rights  reserved. 


*  Corresponding  author.  116  Manning  Dr.,  Mary  Ellen  Jones  Building  906a,  Chapel  Hill,  NC  27599,  USA.  Tel.:  +1  919  845  5585;  fax:  +1  919  962  2388. 
E-mail  address:  ssoper@email.unc.edu  (S.A.  Soper). 

http://dx.doi.org/10.1016/jjpowsour.2014.06.136 

0378-7753 /©  2014  Elsevier  B.V.  All  rights  reserved. 


J.M.  Jackson  et  al.  /  Journal  of  Power  Sources  269  (2014)  274—283 


275 


Nomenclature 

AAP 

change  in  pressure  drop  (Pa) 

R 

resistance  (kg  m~2  s-1) 

A 

cross-sectional  area  of  channel  (m2) 

n 

resistance  of  ith  parallel  channel  (kg  m“2  s_1) 

di 

cross-sectional  area  of  ith  parallel  channel  (m2) 

Ri 

resistance  of  ith  inlet  header  (kg  nrT2  s-1) 

Qp 

cross-sectional  area  of  parallel  channels  (m2) 

R'i 

resistance  of  ith  outlet  header  (kg  m-2  s_1) 

A, 

cross-sectional  area  of  ith  inlet  header  (m2) 

Rj 

resistance  of  tee  junction  (kg  m~2  s_1) 

a; 

cross-sectional  area  of  ith  outlet  header  (m2) 

Re 

Reynolds  number 

An 

cross-sectional  area  of  plate  inlet  (m2) 

U 

velocity  field  (m  s-1) 

Aj 

contact  area  of  tee  junction  (m2) 

Vi 

velocity  of  ith  parallel  channel  (m  s-1) 

D 

hydraulic  diameter  (m) 

Vp 

velocity  of  parallel  channels  (m  s-1) 

F 

friction  factor 

V, 

velocity  of  ith  inlet  header  (m  s-1) 

F\ 

non-uniformity  index 

V'i 

velocity  of  ith  outlet  header  (m  s'1) 

f2 

asymmetry  index 

Vin 

velocity  of  plate  inlet  (m  s_1) 

A 

volumetric  drag  force  (N  m  3) 

w 

width  (m) 

H 

channel  height  (m) 

Wi 

width  of  ith  inlet  header  (m) 

L 

channel  length  (m) 

N 

number  of  channels 

Greek  letters 

P 

channel  perimeter  (m) 

a 

aspect  ratio 

A  P 

pressure  drop  (Pa) 

viscosity  (kg  m-1  s_1) 

1.  Introduction 

Proper  design  of  gas  distributors  for  planar  fuel  cells  is  critical  to 
realize  optimal  operation  and  maximum  power  output  of  a  fuel  cell 
stack.  To  date,  many  flow  designs  have  been  proposed,  evaluated, 
and  are  being  used  for  reagent  delivery  in  commercial  fuel  cells 
[1-4  .  Among  different  configurations,  two  extremes  can  be 
defined:  (i)  A  single  serpentine  channel  covering  the  entire  elec- 
trochemically  active  area  of  the  fuel  cell;  and  (ii)  an  array  of  parallel 
channels,  herein  focusing  on  designs  with  inlet  and  outlet  channels 
arranged  in  the  so-called  Z-type  configuration  [1,5-7  .  The 
serpentine  channel  provides  the  most  uniform  flow  of  reagents,  but 
at  the  same  time,  suffers  from  the  highest  overall  pressure  drop. 
This  is  very  undesirable  for  many  high  power  systems  with  large 
electrode  areas  as  the  pump's  power  parasitically  feeds  from  the 
fuel  cell's  power  output  [5].  On  the  other  hand,  parallel  channel 
configurations  offer  the  lowest  pressure  drop,  as  much  as  an  order 
of  magnitude  lower  than  a  serpentine  channel  covering  the  same 
area  [8],  but  may  suffer  from  severe  flow  maldistribution,  where 
the  middle  channels  are  starved  of  reagent  flow.  This  phenomenon 
can  adversely  affect  a  significant  portion  of  the  electroactive  surface 
area  and  severely  hampers  power  production  of  parallel-configured 
fuel  cells  [2-5,7-20]. 

While  several  methods  have  been  developed  to  numerically 
predict  the  flow  distribution  of  reagents  in  parallel  configurations 


[21  ,  there  are  relatively  few  studies  directed  at  reducing  flow 
maldistribution  by  altering  the  Z-configuration’s  geometry 
[3,18,22-24  .  As  flow  non-uniformity  scales  with  the  number  of 
channels  placed  in  parallel,  various  designs  utilizing  serially 
connected  subsets  of  Z-geometries  with  fewer  parallel  channels, 
termed  discontinuous  geometries,  have  been  proposed.  However, 
applying  a  discontinuous  parallel  configuration  typically  in¬ 
creases  parasitic  pressure  relative  to  a  purely  parallel  design  3], 
and  most  critically,  the  fundamental  non-uniform  profiles  remain 
even  if  they  are  reduced  in  magnitude.  This  last  point  was 
addressed  by  Zhang  et  al.,  who  successfully  corrected  the  non- 
uniform  flow  profiles  in  a  Z-type  configuration  by  adjusting  the 
parallel  channels'  widths  to  increase  flow  through  the  middle 
channels  [18  .  However,  this  optimization  method  was  presented 
in  an  ad  hoc  fashion  as  the  authors  did  not  provide  a  universal 
solution  that  can  be  applied  to  any  geometry,  an  important  point 
considering  parameters  vary  from  one  study  to  the  next.  More¬ 
over,  Kumar,  et  al.  [25]  has  shown  that  optimal  geometric  pa¬ 
rameters  exist  for  the  electroactive  channels  [5].  Altering  the 
widths  of  the  parallel  channels  may  evenly  distribute  reagent 
flow  to  the  middle  channels  but  at  the  expense  of  decreased 
reaction  efficiency.  In  contrast,  header  widths  have  been  adjusted 
to  curb  flow  non-uniformity  in  Z-type  [23 [  and  pin-type  [24] 
configurations,  but  to  date,  no  universal  equation  has  been  pre¬ 
sented  to  predict  how  header  shape  should  be  adjusted  to 
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Fig.  1.  (A)  Schematic  diagram  and  (B)  discrete  representation  of  a  Z-configuration  geometry.  Arrows  represent  the  direction  of  flow. 
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maximize  flow  uniformity  in  any  Z-type  design  without  the  use 
of  an  optimization  algorithm. 

Herein,  we  present  a  simple  and  universal  geometry  optimi¬ 
zation  method  to  ensure  uniform  flow  distribution  in  a  Z-type 
fuel  cell.  We  first  generalize  a  discrete  model  (also  termed  a 
network  analysis  model)  for  Z-type  configurations  [18],  where  the 
geometry  was  subdivided  into  a  network  of  individually  defined 
fluidic  resistors  (see  Fig.  1).  The  simplicity  of  such  a  model  is 
desirable  because  pressure  and  mass  balance  equations  may  be 
defined  and  solved  using  only  linear  algebraic  transformations  of 
matrices,  albeit  these  solutions  are  generally  less  descriptive  than 
analytic  methods  or  CFD  simulations  [18,26  .  We  exploit  this 
simplicity  to  optimize  the  Z-configuration  geometry  by  assuming 
the  case  where  flow  is  perfectly  distributed  throughout  the  fuel 
cell  and  the  geometric  parameters  of  the  parallel  channels  are 
constant.  This  permitted  us  to  reduce  the  governing  pressure  and 
mass  balance  equations  into  straight  forward  geometric  ratios 
between  the  inlet  headers'  resistances  and  areas  that  can  be 
satisfied  by  increasing  header  widths  using  a  simple  algorithm, 
thereby  offering  superior  scalability  to  previously  demonstrated 
optimization  algorithms  [23,24  .  We  demonstrate  this  optimiza¬ 
tion  process  for  several  geometries  and  affirm  the  validity  of 
these  results  by  CFD  simulation.  This  method  comprised  a  simple 
yet  effective  approach  for  elucidating  novel  Z-type  configurations 
of  fuel  cells,  as  well  as  fuel  cell  stacks,  with  markedly  uniform 
flow  distributions. 

As  an  added  benefit,  the  increased  header  widths  reduce 
Reynolds  numbers  throughout  the  headers,  thereby  reducing  flow 
recirculation  in  the  branching  tee  junctions.  In  turn,  this  mini¬ 
mized  asymmetric  non-uniformity  in  reagent  flow,  parasitic  minor 
loss  pressure  drops,  and  reagent  imbalance  between  the  cathode 
and  anode,  all  of  which  have  been  detailed  as  concerns 
[3,22,26,27  .  Furthermore,  we  provide  avenues  by  which  an  opti¬ 
mized  fuel  cell's  footprint  can  be  reduced  while  retaining  flow 
uniformity. 

2.  Computational  methods 

2.2.  Assumptions 

Both  discrete  and  CFD  models  make  the  same  assumptions:  (1) 
The  fluid's  density  (p)  and  viscosity  (/u)  are  constant;  (2)  tempera¬ 
ture  is  held  constant  at  293.15  I<;  (3)  the  fluid  flow  is  steady-state 
and  laminar;  and  (4)  mass  transfer  with  the  electrode  layer  is 
neglected  [18,20]. 


resistance  term  (Rj)  to  a  header  segment  due  to  the  contact  area  of 
the  tee  junction  (At): 


1  (Re  f)pAj 

2  DA 


In  Eq.  (2),  all  geometric  parameters  take  into  account  the 
associated  header  segment  [18  .  Also  note  that  this  model  ne¬ 
glects  any  added  tee  junction  resistance  term  to  the  parallel 
channels  due  to  their  length.  This  assumption  is  addressed  in 
Section  6. 

Using  the  resultant  resistances,  we  can  solve  a  set  of  pressure 
and  mass  balance  equations  to  determine  the  flow  distribution 
through  the  parallel  channels  of  a  Z-configuration  fuel  cell.  The 
pressure  and  mass  balance  equations  for  the  Z-configuration  pre¬ 
sented  in  Fig.  1  are  as  follows: 


v\t\  -T  -T  V2^2 


(3a) 


V2r2  +  ^2^2  =  ^2^2  +  ^3r3 


(3b) 


VinAn  —  V\A  +vlfll 


(3c) 


V\A\  =  V2A2  +  V2a2 


(3d) 


^2^2  =  v3a3 


(3e) 


VftAi  =  V1A1  +  VM  =  v2a2  +  v^A'2  (3f) 

where  v/,  q,  and  a/  are  the  average  linear  velocity,  resistance,  and 
cross-sectional  area  of  the  ith  parallel  channel.  Similarly,  V*,  Rz,  and 
At  describe  the  ith  inlet  header  segment,  and  V/,  R/,  and  A{ 
correspond  to  the  ith  outlet  header  segment. 

Because  each  segment  is  independently  defined,  it  is  rela¬ 
tively  straightforward  to  extend  the  geometry  in  Fig.  1  to  a 
system  with  N  parallel  channels  and  N  -  1  inlet  and  outlet 
header  segments.  Moreover,  we  express  Eq.  (3)  in  matrix  form  so 
that  a  custom  computer  algorithm  can  be  used  to  solve  for  the 
system's  average  linear  velocities  and  flow  distribution.  For 
brevity,  we  apply  both  transformations  simultaneously  and  find 
the  following; 

[M]  =  [R][V]  (4a) 

where 


2.2.  Discrete  model  for  Z-configuration  geometries 

We  adopted  a  discrete  model  [18  in  which  a  Z-type  configu¬ 
ration  fuel  cell  was  segmented  into  a  system  of  individual,  inter¬ 
connected  fluidic  resistors  as  seen  in  Fig.  1.  A  channel  segment's 
geometry  dictates  resistance  (R)  to  fluid  flow  by; 


^in^l  A'n 

^in^An 

A' 


_  1  (Ref)pPL 

2  DA  [  J  [M] 

where  the  channel's  geometry  is  defined  by  its  cross-sectional  area 
(A),  perimeter  (P),  length  (L),  width  (W),  and  height  (H).  The  hy¬ 
draulic  diameter  (D)  is  given  by  4WH/2(  W+H),  and  the  product  of 
the  Reynolds  number  and  friction  factor  (Re  f)  is  approximated  by 
Kays  and  Crawford  [28]  as  13.84+10.38exp(-3.4/a),  where  a  is  the 
channel's  aspect  ratio  (>1).  Additionally,  merging  the  parallel 
channels  and  header  segments  in  a  tee  junction  implicitly  adds  a 


^in^N-lAn 

Avi 

AnVin 

0 


0 


(4b) 
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r\  —r2  0  0 

0  r2  -r3  0 


0 


r'M 
4  ) 


0  0 


0  0 


A'2  J 


0 


0 

0 


0  0 

a!  0  0 

0  a2  0 

0  0  a3 


0  rN-l  ~rN 

0  •••  0 

0  •••  0 

0  ...  0 


0 

A\ 

—A\ 

0 


0  0 
0  0 


aN- l  0  : 

0  aN  0 


0  •••  0 

o  .  o 

A2  : 

~A2  : 

0  0 

:  -An -2  An_  i 

0  ...  0  — A]y_  i 


(4c) 


"1 

v2 

VN 

Vi 

V2 

Vn- i 


(4d) 


Note  that  the  matrices  in  Eq.  (4b  and  c)  reduce  to  those  pub¬ 
lished  by  Zhang  et  al.  [18  if  the  proper  substitutions  are  made. 

We  constructed  an  algorithm  using  the  FORTRAN  77  program¬ 
ming  language  to  solve  Eq.  (4).  Using  any  set  of  fluidic  and  geo¬ 
metric  parameters,  the  program  automatically  assembled  the  [M] 
and  [R]  matrices  using  Eq.  (4b  and  c)  and  inverted  the  [R]  matrix  to 
form  the  [R]_1  matrix.  The  stability  of  this  inversion  was  tested  by 
outputting  the  diagnostic  det([R][R]_1)  =  det([/]),  which  was  unity 
for  systems  where  N  =  1000.  Finally,  [R]_1[M]  was  computed  to 
give  the  [V]  matrix  in  Eq.  (4d  .  The  first  N  elements  of  the  [V]  matrix 
yielded  the  flow  distribution  throughout  their  parallel  channels  of 
the  specified  Z-configuration  geometry,  which  was  characterized 
by  an  established  non-uniformity  index  [8,12] : 

_  max(tg . . .vN)  -  min(iq ...vN) 

1  —  max(v^...vN) 


the  flow  is  perfectly  distributed  through  the  parallel  channels,  i.e., 
vi...vjv  =  vp  in  Eqs.  (3a— e).  We  then  recognize  that  the  pressure 
balance  equations  in  Eq.  (3a  and  b)  simplify  to; 

ViRi  =  Vft  (6a) 

and  the  mass-balance  equations  in  Eqs.  (3c— e)  imply  that; 


\7]yA]y  —  VpQp 

(6b) 

\7jV-iAjv_i  =  VnAn  +  VpQp  =  2  vpQp 

(6c) 

which  can  be  generalized  for  the  ith  inlet  header  by; 


(JV  +  1  — 1> 


CLpVp 


(6d) 


We  now  make  a  second  assumption  that  the  entire  Z-configu¬ 
ration  is  symmetric  (symmetry  plane  shown  in  Fig.  1A).  Here,  the 
geometry  of  the  ith  inlet  header  is  identical  to  the  (N  -  i)th  outlet 
header.  For  the  geometry  in  Fig.  1,  this  means  that  A\  =  A3'  and 
Ax  =A3,  which  is  also  true  for  widths,  heights,  resistances,  and 
average  velocities.  By  Eq.  (6a),  this  constraint  immediately  leads  to 
a  universal  set  of  solutions  by  relating  the  ith  and  (N  +  1  -  i)th  inlet 
headers; 


In  the  case  of  perfectly  uniform  distribution,  ¥\  =  0,  and  F\  -►  1 
as  flow  non-uniformity  becomes  increasingly  severe. 


V^i  =  V'R'i  =  VN+,_tRN+, 


-l 


(6e) 


We  then  substitute  Eq.  (6d)  on  both  sides  of  Eq.  (6e)  to  give; 


2.3.  Discrete  geometry  optimization 

To  optimize  a  Z-type  parallel  configuration  of  a  fuel  cell  for  flow 
uniformity,  we  must  make  a  few  assumptions.  First,  we  assume  that 


/at  .  -i  A  cipVpRi  ■  aPvPRN+i-i 

(n  +  i  -  iy—x — -  *■— * - 

Ai  AN+\-i 

which  reduces  to 


Table  1 

Non-optimized  geometric  parameters  (units:  mm),  F\  parameters  from  discrete  and  CFD  results,  CFD  pressure  drops  (units:  Pa),  and  references  for  geometries  I— III. 


Geometry 

Discrete 

CFD 

Parallel  channel  properties 

Rib  width 

Inlet  width 

Inlet  height 

Ref. 

Fi 

Fi 

A  P 

Number 

Length 

Width 

Height 

I 

0.40 

0.41 

11.1 

11 

50.00 

1.50 

0.60 

1.50 

3.00 

0.60 

[18] 

II 

0.78 

0.79 

30.1 

21 

50.00 

1.50 

0.60 

1.50 

3.00 

0.60 

[18] 

III 

0.92 

0.92 

39.6 

26 

50.00 

2.00 

0.72 

2.00 

4.00 

0.72 

[8] 
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Width  Width 


Fig.  2.  Labeled  drawing  of  geometry  I.  Refer  to  Tables  1  and  3  for  all  parameters  of 
geometries  I — III  and  IV-V,  respectively. 


the  same  dimensions  as  the  parallel  channels,  which  propagates  to 
an  optimized  geometry  with  a  minimal  footprint,  unless  the  N/2... 
N- 1  headers  are  narrowed  further.  It  should  be  noted  that  if  the 
widths  of  the  N/2... N-l  headers  are  increased,  all  widths  will  in¬ 
crease,  and  this  technique  can  be  used  as  a  tool  to  reduce  Reynolds 
numbers  throughout  the  headers  and  the  corresponding  minor 
losses,  distribution  asymmetry,  and  parasitic  pressure  drops  illus¬ 
trated  in  Section  6.  In  other  applications,  we  have  written  algo¬ 
rithms  to  limit  parameters  such  as  Reynolds  number  or  fluidic 
shear  stress  throughout  the  header  channels  by  stipulating  that  Eq. 
(8)  must  be  satisfied  and  flow  through  all  header  satisfies  the  sec¬ 
ondary  condition.  Practically,  this  is  implemented  by  wrapping  the 
width  optimization  algorithm  detailed  later  in  this  section  by  a 
similar  algorithm  stipulating  the  secondary  condition.  These  ap¬ 
plications  were  equally  successful  and  resulted  in  alternative,  more 
linear,  header  shapes  unique  to  these  restrictions. 

Next,  the  geometries  of  the  first  1...N/2  inlet  headers  are  set  to 
satisfy  the  relation  in  Eq.  (8),  which  involves  increasing  their  width 
and/or  height  relative  to  the  last  N/2... N-l  inlet  headers.  Note  that 
if  there  are  an  odd  number  of  parallel  channels,  the  middle  N/2  inlet 
header,  which  is  not  related  to  any  other,  is  simply  assigned  the 
geometry  of  the  adjacent  N/2+1  inlet  header. 

Since  Eq.  (9)  is  far  from  a  simple  expression  relating  W  to  R/A,  it 
is  not  trivial  to  fit  a  universal  expression  to  approximate  the  set  of 
Wi  for  any  geometry  since  there  are  also  dependencies  on  channel 
height  and  length.  Rather,  we  wrote  a  simple  search  algorithm  to 
find  Wi  for  the  first  1...N/2  inlet  headers  by  the  following 
operations: 


Ri  _  *  fijV+l-i  /m 

At  (N  +  1  -  i)  An+1_i  {  } 

The  R/A  ratio  for  a  header  segment  is  given  by  simplification  of 
Eqs.  (1)  and  (2); 


ixL  I  13.84  +  10.38  e-Tr  J  (W  +  H)1 2 
A  =  2  W3H3 

/xAt  M3. 84  +  10.38-e-Tr 

+  4-W3H3 

assuming  that  the  segment's  aspect  ratio  is  given  by  W/H,  i.e.,  W/ 
H>  1. 

Eq.  (8)  is  the  key  to  our  optimization  method  because  it  relates 
the  first  1...N/2  inlet  headers  to  the  last  N/2... N-l  inlet  headers. 
Practically,  we  can  change  the  dimensions  of  the  1...N/2  inlet 
headers  with  respect  to  the  N/2... N-l  inlet  headers  to  satisfy  Eq. 
(8).  It  is  entirely  possible  to  vary  header  heights  to  satisfy  Eq.  (8), 
but  this  approach  could  generate  several  problems.  If  the  parallel 
channels  heights  are  not  altered  to  match  the  header  heights,  re¬ 
agent  must  flow  over  an  additional  surface  from  deeper  inlet 
headers  to  shallower  parallel  channels,  and  this  could  induce  sig¬ 
nificant  minor  losses  at  higher  Reynolds  numbers  similar  to  those 
observed  in  Section  6  of  this  publication  and  unknown  effects  on 
the  flow  distribution.  However,  if  the  parallel  channel  heights  are 
altered,  this  could  negatively  impact  optimal  power  generation  in 
the  parallel  channels  as  shown  by  Kumar  et  al.  5,25].  Therefore,  we 
chose  to  alter  the  headers'  widths  to  satisfy  Eq.  (8). 

There  are  infinite  solutions  to  the  outlined  R/A  relationships; 
some  constraint  must  be  emplaced  to  arrive  at  a  unique  solution, 
which  we  arbitrarily  assigned  as  a  minimized  footprint  herein. 
Thus,  we  begin  by  setting  the  widths  of  the  N/2 . .  .N- 1  headers  with 


(1)  For  the  N/2... N-l  inlet  headers,  calculate  the  R/A  ratios,  and 
set  the  target  R/A  ratios  for  the  1...N/2  inlet  headers  by  Eq. 
(8).  All  subsequent  steps  regard  the  1...N/2  inlet  headers. 

(2)  Set  Wi  equal  to  zero,  and  an  initial  step  size  of  1  mm. 

(3)  Add  a  step  to  the  initial  W\  and  recalculate  R//A;  via  Eq.  (9). 

(4)  If  the  new  R//A,  is  greater  than  the  result  from  step  (1),  add 
another  step  increment  and  repeat  step  (2). 

(5)  If  the  new  RijAi  is  less  than  the  target  from  step  (1 ),  take  one 
step  back  and  decrease  the  step  size  by  a  factor  of  10.  Proceed 
with  step  (2)  unless  the  new  step  size  is  less  than  a  specified 
tolerance  increment.  If  the  tolerance  limit  has  been  reached, 
compare  Rj/A,  as  is  to  Rj/A,  with  an  added  tolerance  step,  and 
choose  the  value  closest  to  the  target.  In  this  study,  we 
specified  the  tolerance  increment  at  0.01  mm  to  reflect 
fabrication  limits  [2  . 

(6)  Set  the  ith  outlet  header  width  equal  to  the  (N+l-i)th  inlet 
header  width. 

The  program  described  in  Section  2.2  was  modified  to  include 
this  search  algorithm  and  solve  for  both  the  initial  and  optimized 
flow  profiles. 

2.4.  CFD  modeling 

As  a  matter  of  validation,  we  used  COMSOL  Multiphysics®  4.3a 
to  conduct  CFD  simulations  of  both  air  and  hydrogen  flow  distri¬ 
butions  through  Z-type  configuration  designs.  Geometries  were 
constructed  within  COMSOL  as  two-dimensional  to  ensure  nu¬ 
merical  tractability.  To  account  for  this  approximation,  a  volumetric 
drag  term  was  added  to  the  velocity  field  (u),  =  -12 juu/H2  .  The 


1  We  analyzed  a  175  channel  Z-configuration  with  an  Intel  i7-3517U  CPU  in  only 

2.625  s  CPU  time.  Moreover,  with  UNC-Chapel  Hill's  KillDevil  supercomputing 

cluster  running  an  Intel  Xeon  X5650  CPU,  we  analyzed  a  1000  channel  configura¬ 
tion  in  376.060  s  of  CPU  time. 


( W  +  H ) 


(9) 
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Fig.  3.  Air  flow  distributions  in  (A)  initial  and  (B)  optimized  geometries  I — III  calculated  via  (green)  discrete  and  (red)  CFD  methods.  Dashed  lines  represent  perfectly  distributed 
profiles.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 
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Fig.  4.  CFD  models  of  air  flow  in  (A)  initial  and  (B)  optimized  geometries  I — III. 
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validity  of  this  volumetric  drag  term  was  affirmed  by  comparison 
with  a  three-dimensional  model  of  geometry  I  defined  in  Section  3 
(data  not  shown).  Elongated  inlets  and  outlets  were  used  to  sta¬ 
bilize  flow  profiles  prior  to  flow  encountering  the  parallel  channels. 

A  faux  internal  boundary  was  drawn  across  the  middle  of  the 
parallel  channels.  This  boundary  had  no  effect  except  to  ensure  that 
the  meshing  algorithm  assigned  data  points  along  this  boundary  in 
each  channel,  which  were  used  to  construct  velocity  line  plots. 
Solutions  were  obtained  via  meshing2  and  solving3  with  custom 
settings  within  COMSOL.  Excluding  points  that  defined  the  wall, 
data  was  averaged  to  generate  the  linear  velocity  through  a  parallel 
channel.  To  account  for  small  deviations  caused  by  extracting  the  vz- 
data  in  this  manner,  sets  of  vz  from  both  discrete  and  CFD  solutions 
were  normalized  with  respect  to  the  average  linear  velocity  over  all 
V/.  Pressure  drops  were  calculated  using  two  lines  across  the  inlet 
and  outlet  channels,  which  were  directly  adjacent  to  the  first  and 
last  parallel  channel,  respectively,  to  account  for  the  elongated  in¬ 
lets  and  outlets.  The  maximum  pressure  of  the  inlet  line  was  sub¬ 
tracted  from  the  minimum  of  the  outlet  line. 

Optimized  geometries  from  discrete  solutions  in  Section  2.3 
were  constructed  within  COMSOL  by  fitting  cubic  functions 
through  each  header  segment  at  the  branching  tee  junctions.  This 
resulted  in  a  smooth  transition  between  the  widths,  thereby 
avoiding  abrupt  changes  in  fluid  flow  and  minor  losses  due  to 
sudden  contractions.  It  must  be  noted  that  while  the  cubic  function 
resulted  in  well  distributed  flow  (see  Section  4),  this  adaptation  was 
empirical.  It  is  entirely  plausible  that  there  are  alternative  methods 
of  fitting  the  sets  of  header  widths  that  better  match  the  discrete 
optimization  results.  Moreover,  potential  deviations  in  the  fabri¬ 
cation  of  these  curvatures  and  their  impact  on  flow  distribution 
warrant  future  experimental  validation. 

3.  Validation 

Using  both  discrete  (Section  2.2)  and  CFD  (Section  2.4)  methods, 
we  assessed  air  flow  distributions  for  three  published  Z-configu- 
ration  geometries  [8,18]  that  are  parameterized  in  Table  1  and  Fig.  2. 
For  comparability,  all  inlet  velocities  were  set  so  that  the  average 
parallel  channel  velocity  should  be  0.1  m  s-1  [8]  if  perfectly 
distributed,  i.e.,  Vin  =  Eiliai) /An -0.1  ms-1  (stoichiometry  ratio  of 
1.3  and  3.3  for  air  and  hydrogen,  respectively,  at  a  current  density  of 
0.3  A  cm-2). 

Flow  distributions  of  geometries  I — III  are  shown  in  Fig.  3A,  and 
CFD  velocity  surfaces  are  shown  in  Fig.  4A.  Both  sets  of  F\  param¬ 
eters,  calculated  by  Eq.  (5)  and  shown  in  Table  1,  are  nearly  identical 
between  discrete  and  CFD  results.  Additionally,  these  results  coin¬ 
cide  with  previously  published  results  [8,18  ,  all  of  which  validates 
the  discrete  method  for  these  geometries. 

4.  Geometry  optimization 

We  applied  the  discrete  optimization  code  in  Section  2.3  to 
geometries  I— III  and  adapted  the  resultant  inlet  and  outlet  header 
widths  to  the  CFD  simulations  as  described  in  Section  2.4.  The  F\ 
parameters,  inlet  widths,  and  percent  changes  in  the  pressure 


2  The  maximum  element  size,  minimum  element  size,  and  maximum  element 
growth  rate  were  0.25  mm,  10.3  pm,  and  1.04,  respectively;  the  resolutions  of 
curvature  and  narrow  regions  were  0.1  and  16,  respectively.  The  geometries  pre¬ 
sented  herein  consisted  of  approximately  150,000  to  600,000  elements. 

3  Systems  were  solved  using  the  PARDISO  algorithm,  the  Double-Dogleg 
nonlinear  solver,  automatic  pseudo-time  stepping,  and  a  relative  tolerance  that 
was  minimized  for  each  geometry  to  ensure  convergence  to  a  unique  solution  in  all 
cases.  For  the  largest  25-channel  geometry,  it  took  4111  s  of  CPU  time  using  an  Intel 
i7-3770K  processor. 


required  to  drive  the  system  from  CFD  solutions  are  shown  in 
Table  2.  The  optimized  flow  profiles  from  discrete  and  CFD  solu¬ 
tions  are  shown  in  Fig.  3B,  and  the  CFD  velocity  profiles  are  pre¬ 
sented  in  Fig.  4B. 

After  optimization,  we  have  significantly  reduced  flow  maldis¬ 
tribution  in  geometries  I— III,  where  the  air  flow  F\  parameters 
decreased  by  86%  on  average.  Additionally,  the  parasitic  pressure 
required  to  drive  these  optimized  geometries  is  either  slightly 
reduced  or  essentially  unaffected  (see  Table  2).  Thus,  this  optimi¬ 
zation  method  is  effective,  simply  devised  for  any  given  system  via 
the  relations  in  Section  2.3,  does  not  change  any  geometric  pa¬ 
rameters  of  the  parallel  channels  that  could  affect  reaction  effi¬ 
ciencies,  and  offers  improved  scalability  over  previous  optimization 
methods  [23,24],  which  is  immediately  evident  from  the  CPU  times 
in  Sections  2.3  and  2.4  that  were  required  to  solve  for  velocity  fields 
and  would  be  further  improved  by  more  advanced  algorithms  than 
the  simplified  search  algorithm  outlined  in  Section  2.3. 

We  recognize  that  the  single-phase,  steady-state  models 
employed  for  geometry  optimization  are  limited  in  scope  and  do 
not  describe  dynamic  processes,  such  as  the  generation  of  water 
droplets/slugs  and  water  purging  via  hydrodynamic  force,  that 
create  time-dependent  fluctuations  in  pressure  and  flow  distribu¬ 
tions  during  device  operation.  Fuel  cell  flooding  remains  a  critical 
limitation  of  parallel  systems,  such  as  Z-type  geometries,  primarily 
since  flow  velocities  and  hydrodynamic  forces  are  significantly 
reduced  (relative  to  single,  serpentine  channels)  when  stoichiom¬ 
etry  is  conserved  [29-32  .  The  uniform  steady-state  flow  fields 
illustrated  in  this  section  are  also  necessary  to  provide  uniform 
hydrodynamic  forces  for  water  removal  across  the  entire  electro¬ 
active  surface.  As  such,  the  model  presented  herein  functions  as  a 
necessary  initial  step  in  the  development  of  more  complex  multi¬ 
phase  models  describing  water  removal  dynamics. 

5.  Minimizing  footprint  of  optimized  geometries 

As  can  be  seen  in  Fig.  4,  geometry  Ill's  footprint  increased  sub¬ 
stantially  after  optimization.  While  the  optimized  header  width  is 
on  the  order  of  stack  manifold  headers  [33]  and  likely  not  an  issue 
for  fabrication,  it  is  inevitable  that  during  scale-up,  the  optimized  Z- 
configuration’s  footprint  could  increase  beyond  tolerable  limits  or 
the  stack  manifolds  could  become  too  large.  To  address  this  po¬ 
tential  issue,  we  evaluated  several  avenues  by  which  the  leading 
header  widths  and  total  footprint  can  be  minimized  without  per¬ 
turbing  the  uniform  flow  in  the  parallel  channels. 

First,  rather  than  applying  the  exact  widths  satisfying  Eq.  (8),  we 
reduced  the  headers  of  the  optimized  geometry  III  (Fig.  4B)  by 
different  percentages  of  the  optimal  widths  to  determine  if  we 
could  apply  a  less  exact  geometry  optimization  with  a  smaller 
footprint  but  still  yield  a  uniform  flow  distribution.  The  discrete 
results  from  applying  various  percentages  of  the  optimal  widths  to 
geometry  III  are  shown  in  Fig.  5.  Note  that  these  percentages  regard 
the  increases  in  widths  of  the  1 . .  .N/2  headers  relative  to  the  widths 
of  the  N/2...N-1  headers  (see  Section  2.3).  Applying  80%  of  the 

Table  2 

Optimized  geometries'  inlet  widths  (units:  mm),  parameters  calculated  via 
discrete  and  CFD  analyses,  and  CFD  pressure  drops  (units:  Pa)  and  their  percent 
decrease  relative  to  initial  geometries. 


Geometry 

Inlet  width 

Discrete 

CFD 

Fi 

Fi 

A  P 

AAP  (%) 

I 

12.28 

0.00 

0.06 

10.4 

-6.8 

II 

24.37 

0.00 

0.09 

30.0 

-0.3 

III 

41.72 

0.01 

0.14 

39.3 

-0.6 
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Fig.  5.  Discrete  Fi  parameters  for  optimizations  of  geometry  III  using  various  per¬ 
centages  of  the  optimal  header  widths.  The  line  connecting  data  points  is  for  illus¬ 
tration  only. 


fundamental  maldistribution  profiles  are  still  present.  Even  smaller 
discontinuous  subsets  of  geometry  III  would  be  necessary  to  further 
reduce  the  magnitude  of  non-uniformity,  but  again,  this  optimi¬ 
zation  technique  would  result  in  increased  parasitic  pressure  [3]. 
However,  if  we  connect  optimized  subsets  of  geometry  III  in  a 
discontinuous  fashion  (Fig.  6B),  maldistribution,  pressure,  and 
overall  footprint  is  reduced. 

Lastly,  channel  heights  could  be  increased  throughout  the  entire 
geometry  to  minimize  the  optimized  geometry's  overall  footprint. 
This  necessitates  smaller  increases  in  header  widths  to  produce 
larger  decreases  in  the  R/A  ratio  in  Eq.  (9).  If  geometry  Ill's  depth  is 
increased  to  1.5  mm,  which  has  been  suggested  as  optimal  for 
hydrogen  consumption  [25  ,  the  optimized  leading  inlet  width 
then  decreases  from  41.72  mm  to  30.71  mm.  This  technique  is 
limited  because  further  increases  in  channel  height  could  reduce 
reaction  efficiency  in  the  parallel  channels.  Also,  changes  in  the 
headers'  hydraulic  diameter  and  inlet  linear  velocity  may  increase 
Reynolds  numbers  throughout  the  inlets,  which  can  cause  flow  to 
recirculate  at  the  branching  tee  junctions  and  asymmetry  in  the 
flow  distribution  due  to  the  minor  losses  described  in  Section  6. 


optimal  widths  increased  the  Fi  parameter  from  0.01  to  0.21,  but 
applying  90%  resulted  in  an  F\  of  0.11,  which  is  within  the  ±5% 
deviation  that  has  been  previously  defined  as  an  acceptable  toler¬ 
ance  for  flow  non-uniformity  [22  . 

Second,  discontinuous  designs  may  be  applied  to  optimized 
geometries  in  order  to  reduce  the  fuel  cell's  footprint.  For  example, 
a  discontinuous  geometry  III  is  shown  in  Fig.  6A.  It  is  clear  that  the 


Fig.  6.  CFD  air  flow  profiles  through  geometry  III  in  a  (A)  discontinuous  configuration 
and  an  (B)  optimized  discontinuous  configuration.  The  scale  bar  is  shown  in  Fig.  3. 


6.  Minor  losses 

As  mentioned  in  Section  2.2,  resistances  within  the  discrete 
model  only  concern  major  pressure  drops  to  viscous  drag.  In  cases 
where  Reynolds  numbers  are  large  enough  to  induce  flow  recir¬ 
culation  in  the  branching  tee  junctions,  the  minor  loss  pressure 
drops,  asymmetric  skew  in  flow  distributions,  and  reagent  imbal¬ 
ance  between  the  cathode  and  anode  will  not  be  predicted  by  this 
model  [3,22,26,27  .  We  must  explicitly  state  that  it  is  not  a  trivial 
task  to  describe  minor  losses  occurring  in  branching  and  combining 
tee  junctions  in  the  algorithms  in  Section  2.2.  Minor  effects  are  a 
function  of  the  ratio  of  velocities  between  the  branched  and  com¬ 
bined  flow  [26  ,  whereas  velocities  are  individually  defined  in  the 
[V]  matrix.  Including  minor  effects  would  require  reformulating  the 
entire  discrete  model. 

To  demonstrate  this  issue,  we  evaluated  geometries  IV  and  V 
that  are  parameterized  in  Table  3.  Geometry  IV  was  a  Z-configu- 
ration  adapted  from  channel  dimensions  optimized  for  hydrogen 
consumption  by  Kumar  et  al.  [25  ,  and  geometry  V  was  adapted 
from  a  commercially  available  Z-configuration  fuel  cell  studied  by 
Iranzo  et  al.  [2].  To  achieve  0.1  m  s_1  air  flow  in  the  parallel  chan¬ 
nels  of  geometries  IV  and  V,  the  Reynolds  numbers  at  the  inlets 
were  129  and  140,  respectively,  and  in  Fig.  7,  we  show  flow  recir¬ 
culation  developing  along  the  branching  tee  junctions  (experi¬ 
mentally  demonstrated  by  Barreras  et  al.  19]).  These  results  are 
contrasted  to  geometry  I  in  Figs.  4A  and  7,  where  the  Reynolds 
number  at  the  inlet  was  only  37,  and  no  recirculation  was  observed. 

As  the  fluid's  velocity  and  Reynolds  number  decreases  along  the 
inlet  headers,  recirculation  and  the  accompanying  minor  loss 
pressure  drops  decrease  in  magnitude  [26].  Thus,  parallel  channels 
farther  from  the  inlet  are  biased  with  a  lower  resistance  and  faster 
flow  through  the  parallel  channels,  thereby  modulating  the  discrete 
model's  parabolic  shape  to  generate  the  asymmetric  distributions 
shown  in  Fig.  8A.  Thus,  for  the  non-optimized  geometry  I,  no  sig¬ 
nificant  air  flow  distribution  asymmetry  was  observed,  and  the 


Table  3 

Non-optimized  geometric  parameters  (units:  mm)  and  references  for  geometries 
IV-V. 


Geometry 

Parallel  channel  properties 

Rib 

Inlet 

Inlet 

Ref. 

Number 

Length 

Width 

Height 

width 

width 

height 

IV 

13 

37.00 

1.50 

1.50 

0.50 

1.50 

1.50 

[25] 

V 

25 

72.25 

2.00 

0.80 

0.80 

2.00 

0.80 

[2] 
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Fig.  7.  Inlet  header  velocity  profiles  showing  recirculation  and  minor  losses  in  geometries  I,  IV,  and  V.  Arrows  indicate  direction  of  flow. 
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Fig.  8.  The  (A)  initial  and  (B)  optimized  CFD  flow  distributions  of  (red)  air  and  (blue)  hydrogen  and  (green)  discrete  solution  through  geometries  IV  and  V.  Dashed  lines  represent 
perfectly  distributed  profiles.  (For  interpretation  of  the  references  to  color  in  this  figure  legend,  the  reader  is  referred  to  the  web  version  of  this  article.) 


non-optimized  geometries  II  and  III,  which  have  inlet  Reynolds 
numbers  of  70  and  105,  respectively,  showed  slight  asymmetry 
forming  in  the  air  flow  distributions  (see  Fig.  3A). 

To  quantitate  these  effects,  we  define  a  new  flow  asymmetry 
parameter: 


\ 

\ 

max 

Vi  ...Vn 

2 

-  max 

Vn.  .  .Vi\j 

2 

/ 

/ 

max(iq  ...vN) 


(10) 


In  all  discrete  calculations,  ¥2  will  be  zero  as  the  flow  distribu¬ 
tion  is  always  symmetric.  But  in  CFD  simulations,  ¥2  -►  1  as 
asymmetry  and  minor  losses  become  more  problematic.  For  ge¬ 
ometries  IV  and  V,  both  ¥\  and  ¥2  parameters  are  shown  in  Table  4 
along  with  CFD  pressure  drops.  As  a  note,  geometries  I,  II,  and  III  in 
Fig.  4A  have  air  flow  ¥2  parameters  of  0.03,  0.07,  and  0.10, 
respectively. 

Hydrogen  flow  profiles  are  less  skewed  and  more  closely  match 
discrete  results  since  hydrogen's  kinematic  viscosity  and  Reynolds 


numbers  are  an  order  of  magnitude  less  than  air's.  As  noted  pre¬ 
viously,  this  creates  reagent  imbalance  between  the  cathode  and 
anode  that  can  negatively  impact  a  fuel  cell's  efficiency  [3].  As  ex¬ 
pected,  hydrogen  flow  ¥2  parameters  for  these  geometries  are 
<0.01  for  all  cases. 

Even  though  these  minor  losses  are  not  described  in  the  discrete 
model  detailed  in  Section  2.2  or  the  optimization  method  in  Section 
2.3,  the  increased  header  widths  in  optimized  geometries  inad¬ 
vertently  reduce  these  minor  losses  by  reducing  Reynolds  numbers 
due  to  slower  reagent  flow,  which  outweighs  increasing  hydraulic 
diameters.  For  example,  after  optimizing  geometries  IV  and  V, 


Table  4 

Non-optimized  Fi  parameters  from  discrete  and  CFD  analyses,  and  CFD  air  and 
hydrogen  flow  distribution  F2  parameters  and  pressure  drops  (units:  Pa). 


Geometry 

Discrete 

CFD  (H2) 

CFD  (air) 

Fi 

Fi 

f2 

A  P 

Fi 

f2 

A  P 

IV 

0.67 

0.70 

0.10 

2.2 

0.80 

0.59 

5.7 

V 

0.92 

0.92 

0.03 

21.2 

0.93 

0.31 

48.0 

J.M.  Jackson  et  al.  /  Journal  of  Power  Sources  269  (2014)  274—283 


283 


Table  5 

Optimized  geometries'  inlet  widths  (units:  mm),  Fi  parameters  calculated  via  discrete  and  CFD  analyses,  and  CFD  A P  pressure  drops  (units:  Pa)  and  their  percent  decrease 
relative  to  initial  geometries. 


Geometry 

Inlet  width 

Discrete 

CFD  (H2) 

CFD  (air) 

Ei 

Ei 

E2 

A  P 

AAP  (%) 

Ei 

e2 

A  P 

AAP  (%) 

IV 

9.86 

0.00 

0.06 

0.02 

1.2 

-46 

0.27 

0.21 

2.7 

-53 

V 

41.23 

0.00 

0.10 

0.01 

11.2 

-47 

0.14 

0.05 

23.7 

-51 

pressures  were  reduced  to  about  half  as  the  minor  loss  pressure 
drops  lessened  ( Table  5).  The  F2  parameters  of  geometries  I-III  were 
all  reduced  to  0.02  after  optimization,  and  these  asymmetry  pa¬ 
rameters  of  geometries  IV  and  V  were  also  reduced,  albeit  not  to  the 
same  extent.  As  mentioned  in  Section  2.3,  rather  than  constrain 
optimization  results  by  minimized  footprint  (by  restraining  the  N/ 
2...N-1  headers  to  the  parallel  channel  width),  it  is  quite  straight 
forward  to  sequentially  increase  the  N/2...N-1  headers  until  flow 
through  the  optimized  corresponding  header  has  a  Reynolds 
number  less  than  a  specified  value,  such  as  100. 


7.  Conclusions 

Uniform  delivery  of  reagents  to  the  Z-configuration’s  parallel 
channels  is  imperative  for  optimum  performance  of  a  fuel  cell  stack 
[2-5,7-20  .  We  modified  a  simple  discrete  method  to  assess  flow 
maldistribution  and  deduced  a  mathematical  relationship  for 
increasing  header  widths  to  optimize  the  flow  distribution.  We 
have  shown  the  accuracy  of  the  discrete  method,  the  success  of 
optimization  at  reducing  flow  maldistribution,  and  the  reduction  of 
parasitic  minor  loss  pressure  drops  and  the  corresponding  flow 
distribution  asymmetry.  In  all  of  these  optimizations,  the  parallel 
channel's  geometry,  which  is  often  already  optimized  for  reaction 
efficiency,  is  not  altered,  and  solutions  are  obtained  with  compu¬ 
tational  ease.  To  the  best  of  our  knowledge,  this  represents  the  first 
time  a  model  has  been  described  to  curb  flow  maldistribution  that 
is  universally  applicable  to  all  Z-type  configurations  of  planar  fuel 
cells  [22  .  Furthermore,  this  research  is  equally  significant  in  the 
design  of  fuel  cell  stack  manifolds,  where  reagent  delivery  between 
layers  of  the  stack  exhibits  similar  maldistribution  [34  . 

Water  management  was  not  addressed  in  this  manuscript  and 
remains  a  critical  issue  for  parallel  PEM  fuel  cells.  The  multi-phase 
models  [31,35]  and  empirical  measurements  [29,36]  that  are 
required  to  address  water  management  in  these  optimized  geom¬ 
etries,  especially  with  respect  to  the  water  dynamics  occurring  in 
the  optimized  header  channels,  are  currently  under  investigation. 
In  alternative  fields  of  study,  such  as  microfluidics,  where  large 
scale  Z-configuration  networks  are  needed  for  liquid  delivery 
[37-39],  this  optimization  method  is  immediately  applicable. 
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